require(Ternary)
## Loading required package: Ternary
## Warning: package 'Ternary' was built under R version 4.3.2
Pcol <- rgb(123/255,188/255,73/255)
Vcol <- rgb(168/255,68/255,151/255)
Mcol <- rgb(231/255,186/255,96/255)
# Parameters
# Allocations
alphaV <- 0.3
alphaP <- 0.1
alphaM <- 1-alphaV-alphaP
# External resource availability
L <- 100
I <- 10
B <- 1e6
# Respiratory costs
r <- 0.1 #0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Set up the simulation
tset <- seq(from = 0, to = 10, by = 0.0001)
Vset <- NaN*tset; Vset[1] <- alphaV
Pset <- Vset; Pset[1] <- alphaP
Mset <- Vset; Mset[1] <- alphaM
growth <- NaN*tset
# Run the simulation
for(i in 2:length(tset)){
dt <- tset[i] - tset[i-1]
jC <- dt*( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
jN <- dt*( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
jG <- dt*( yGM*Mset[i-1] )
growth[i] <- min(qC*jC,qN*jN,qG*jG)
Pset[i] <- Pset[i-1] + growth[i]*alphaP
Vset[i] <- Vset[i-1] + growth[i]*alphaV
Mset[i] <- Mset[i-1] + growth[i]*alphaM
}
plot(tset,Pset,type='l',las=1,xlab='Time',ylab='Organelle Structures',ylim=c(0,max(max(Pset),max(Mset),max(Vset))),col=Pcol,lwd=2)
lines(tset,Vset,lwd=2,col=Vcol)
lines(tset,Mset,lwd=2,col=Mcol)
plot(tset,alphaP*growth/Pset,type='l',las=1,col=Pcol,lwd=2,ylim=c(0,.15))
lines(tset,alphaM*growth/Mset,col=Mcol,lwd=2)
lines(tset,alphaV*growth/Vset,col=Vcol,lwd=2)
# Growth rates converge
lm(log(tail(Mset,50))~tail(tset,50))
##
## Call:
## lm(formula = log(tail(Mset, 50)) ~ tail(tset, 50))
##
## Coefficients:
## (Intercept) tail(tset, 50)
## -0.5108 0.6000
lm(log(tail(Vset,50))~tail(tset,50))
##
## Call:
## lm(formula = log(tail(Vset, 50)) ~ tail(tset, 50))
##
## Coefficients:
## (Intercept) tail(tset, 50)
## -1.204 0.600
lm(log(tail(Pset,50))~tail(tset,50))
##
## Call:
## lm(formula = log(tail(Pset, 50)) ~ tail(tset, 50))
##
## Coefficients:
## (Intercept) tail(tset, 50)
## -2.303 0.600
# Ratios converge on allocations
tail(Pset/(Pset+Mset+Vset))
## [1] 0.1 0.1 0.1 0.1 0.1 0.1
tail(Mset/(Pset+Mset+Vset))
## [1] 0.6 0.6 0.6 0.6 0.6 0.6
# Can we predict growth rate?
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
guess_growth <- min(qC*jC,qN*jN,qG*jG)
guess_growth
## [1] 0.6
# It works!!
Show an example of convergence.
# Parameters
# Allocations
alphaV <- 0.3
alphaP <- 0.1
alphaM <- 1-alphaV-alphaP
# External resource availability
L <- 100
I <- 10
B <- 1e6
# Respiratory costs
r <- 0.1#0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Set up the simulation
tset <- seq(from = 0, to = 10, by = 0.0001)
Vset <- NaN*tset; Vset[1] <- 8
Pset <- Vset; Pset[1] <- 3
Mset <- Vset; Mset[1] <- 3.5
growth <- NaN*tset
# Run the simulation
for(i in 2:length(tset)){
dt <- tset[i] - tset[i-1]
jC <- dt*( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
jN <- dt*( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
jG <- dt*( yGM*Mset[i-1] )
growth[i] <- min(qC*jC,qN*jN,qG*jG)
Pset[i] <- Pset[i-1] + growth[i]*alphaP
Vset[i] <- Vset[i-1] + growth[i]*alphaV
Mset[i] <- Mset[i-1] + growth[i]*alphaM
}
Tset <- Pset+Vset+Mset
plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='')
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
Make a growth array
alphaPset <- seq(from=0, to = 1,by = 0.01)
alphaVset <- seq(from = 0, to = 1, by = 0.005)
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
image(growths,xlab='alpha P', ylab='alpha V',las=1)
Make a few showing different outcomes.
alphaPset <- seq(from=0, to = 1,by = 0.01)
alphaVset <- seq(from = 0, to = 1, by = 0.05)
par(mar=c(4,4,1,1),mfrow=c(1,3))
# Option 1: Phototrophy
# External resource availability
L <- 100
I <- 5
B <- 0
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
image(growths,xlab='alpha P', ylab='alpha V',las=1)
# Option 2: Phagotrophy
# External resource availability
L <- 0
I <- 0
B <- 1e6
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
image(growths,xlab='alpha P', ylab='alpha V',las=1)
# Option 3: Mixotrophy
# External resource availability
L <- 100
I <- 5
B <- 1e6
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
#image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(30,'YlOrRd',rev=T))
image(growths,xlab='alpha P', ylab='alpha V',las=1)
OK so looking at just the mixotroph plot we can see that there are three “ridges” to the surface.
par(mar=c(4,4,1,1))
alphaVset <- seq(from = 0, to = 1, length.out = 500)
alphaPset <- seq(from = 0, to = 1, length.out = 495)
# Option 3: Mixotrophy
# External resource availability
L <- 100
I <- 5
B <- 1e6
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
Let’s see if we can trace these.
par(mar=c(4,4,1,1),mfrow=c(1,2))
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
alphaPset.1 <- (qN*uNB*B*alphaVset-qC*(uCB*B*alphaVset-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset)+qC*r-qC*uCB*B*alphaVset)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset) - qN*uNB*B*alphaVset)/(qN*uNI*I + qG*yGM)
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
lines(alphaPset.2,alphaVset,lty=2)
lines(alphaPset.1,alphaVset)
lines(alphaPset.3,alphaVset,lty=3)
legend(x='topright',inset=0.02,legend=c('jC=jN','jC=jG','jN=jG'),lty=c(1,2,3))
Finding the intersection…
plot(alphaVset,abs(alphaPset.1-alphaPset.2),type='l',xlab='alpha V',ylab='delta alpha P')
lines(alphaVset,abs(alphaPset.2-alphaPset.3),col='dodgerblue')
lines(alphaVset,abs(alphaPset.1-alphaPset.3),col='darkgreen')
points(alphaVset[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],0)
#alphaPset.1[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))]
L <- 10
I <- 0
B <- 1e6
growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))
for(i in 1:length(alphaPset)){
for(j in 1:length(alphaVset)){
if(alphaPset[i] + alphaVset[j] <=1){
alphaP <- alphaPset[i]
alphaV <- alphaVset[j]
alphaM <- 1-alphaP-alphaV
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM
growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
}
}
}
par(mar=c(4,4,1,1),mfrow=c(1,2))
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
alphaVset2 <- seq(from = 0, to = 1, length.out = 500)
alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
lines(alphaPset.2,alphaVset2,lty=2)
lines(alphaPset.1,alphaVset2)
lines(alphaPset.3,alphaVset2,lty=3)
points(alphaPset.1[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],alphaVset2[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],pch=21,lwd=1,cex=1.5)
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
points(alphaPstar,alphaVstar,pch=21,col='red',cex=2.25,lwd=1)
legend(x='topright',inset=0.02,legend=c('gC=gN','gC=gG','gN=gG','g* numer.','g* analyt.'),lty=c(1,2,3,0,0),pch=1,pt.cex=c(0,0,0,1.5,2.25),col=c('black','black','black','black','red'))
alphaVstar
## [1] 0
alphaPstar
## [1] 1
alphaMstar <- 1 - alphaVstar - alphaPstar
alphaMstar
## [1] 0
dim(growths)
## [1] 495 500
gstar <- max(growths,na.rm = T)
which(growths==gstar)
## [1] 190081
which(growths==max(growths,na.rm = T), arr.ind=T)
## row col
## [1,] 1 385
alphaVset[which(growths==max(growths,na.rm = T), arr.ind=T)[2]]
## [1] 0.7695391
alphaPset[which(growths==max(growths,na.rm = T), arr.ind=T)[1]]
## [1] 0
(alphaVstar-alphaVset[which(growths==max(growths,na.rm = T), arr.ind=T)[2]])/(alphaVstar)*100
## [1] Inf
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')
# External resource availability
L <- 100
I <- 5
B <- 0
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
Now that we can make these plots, can we add lines to them?
par(mar=rep(0.2, 4))
resolutionsetting <- 100
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)
points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5)
points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5,lty=2)
points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5,lty=3)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),pch=21,col='red',cex=1,lwd=1)
How about a strict strategy?
par(mar=rep(0.2, 4))
resolutionsetting <- 100
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){
alphaVstar <- alphaVstar.1
alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
} else {
alphaVstar <- alphaVstar.2
alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
}
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaPset <- 1-alphaMset-alphaVset2
points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5)
TernaryPoints(c(1-alphaMstar-alphaVstar,alphaMstar,alphaVstar),pch=21,col='red',cex=1,lwd=1)
And now let’s try marking all three.
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')
# External resource availability
L <- 100
I <- 5
B <- 0
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){
alphaPstar <- alphaPstar.1
alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaPstar <- alphaPstar.2
alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
linetypesetting <- 3
}
alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset <- 1-alphaMset-alphaPset2
points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){
alphaVstar <- alphaVstar.1
alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaVstar <- alphaVstar.2
alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
linetypesetting <- 3
}
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaPset <- 1-alphaMset-alphaVset2
points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)
points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=2)
points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty=2)
points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty=3)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
And all three in grayscale.
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')
# External resource availability
L <- 100
I <- 5
B <- 0
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
aP1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
aP2 <- (qG*yGM)/(qN*uNI*I+qG*yGM)
alphaPstar <- max(aP1,aP2)
alphaVstar <- 0
alphaMstar <- 1 - alphaPstar
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
aV1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
aV2 <- (qG*yGM)/(qN*uNB*B+qG*yGM)
alphaVstar <- max(aV1,aV2)
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')
# External resource availability
L <- 100
I <- 5
B <- 0
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){
alphaPstar <- alphaPstar.1
alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaPstar <- alphaPstar.2
alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
linetypesetting <- 3
}
alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset <- 1-alphaMset-alphaPset2
points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){
alphaVstar <- alphaVstar.1
alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaVstar <- alphaVstar.2
alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
linetypesetting <- 3
}
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaPset <- 1-alphaMset-alphaVset2
points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)
points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1)
points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=2)
points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=3)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')
# External resource availability
L <- 100
I <- 5
B <- 0
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){
alphaPstar <- alphaPstar.1
alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaPstar <- alphaPstar.2
alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
linetypesetting <- 3
}
alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset <- 1-alphaMset-alphaPset2
points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')
# External resource availability
L <- 0
I <- 0
B <- 0.88e6 #1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){
alphaVstar <- alphaVstar.1
alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
linetypesetting <- 2
} else {
alphaVstar <- alphaVstar.2
alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
linetypesetting <- 3
}
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaPset <- 1-alphaMset-alphaVset2
points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')
# External resource availability
L <- 10
I <- 2.8 #3 #2
B <- 1e6
# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar
# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)
alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)
points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1)
points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=2)
points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]
TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=3)
TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)
First, let’s make sure we have a robust algorithm for choosing the optimal investment strategy. We’ll confirm using our ternary plot friend.
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,1))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1 # yield from growth
uCL <- 0.01 # carbon from light
uNI <- 1 # nutrients from inorganic nutrients
uCB <- 0.000003 # carbon from bacteria
uNB <- uCB/10 # nitrogen from bacteria
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# External resource availability
L <- 1000
I <- 0
B <- 0.88e6 #1e6
L <- 10
I <- 2.8 #3 #2
B <- 1e6
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
#text(0.,0.95,'Phagotrophy Optimal',col='black')
# Compute the strategy
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
#ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Add analytical solution
strats <- matrix(data=NA,nrow=3,ncol=4)
## Strict phagotrophy
alpha_til_V <- max((qG*yGM+qC*r)/(qC*uCB*B + qG*yGM),(qG*yGM)/(qN*uNB*B+qG*yGM))
alpha_til_P <- 0
alpha_til_M <- 1 - alpha_til_V - alpha_til_P
g_til_CV <- qC*(uCB*B*alpha_til_V-r)
g_til_NV <- qN*uNB*B*alpha_til_V
g_til = qG*yGM*alpha_til_M # growth rate of strict phagotroph
strats[1,1] <- alpha_til_V
strats[1,2] <- alpha_til_P
strats[1,3] <- alpha_til_M
strats[1,4] <- g_til
## Strict phototrophy
alpha_bar_V = 0
alpha_bar_P = max(c((qG*yGM+qC*r)/(qC*uCL*L+qG*yGM),qG*yGM/(qN*uNI*I+qG*yGM)))
alpha_bar_M = 1-alpha_bar_P
g_bar_CP = qC*(uCL*L*alpha_bar_P - r)
g_bar_NP = qN*uNI*I*alpha_bar_P
g_bar = qG*yGM*alpha_bar_M # growth rate of strict phototroph
strats[2,1] <- alpha_bar_V
strats[2,2] <- alpha_bar_P
strats[2,3] <- alpha_bar_M
strats[2,4] <- g_bar
## Mixotrophy
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alpha_star_V <- (term.c-term.d)/(term.a+term.b)
alpha_star_P <- (qG*yGM*(1-alpha_star_V)+qC*r-qC*uCB*B*alpha_star_V)/(qG*yGM+qC*uCL*L)
alpha_star_M <- 1 - alpha_star_V - alpha_star_P
g_star = qG*yGM*alpha_star_M
strats[3,1] <- alpha_star_V
strats[3,2] <- alpha_star_P
strats[3,3] <- alpha_star_M
strats[3,4] <- g_star
## Picking optimum
strats_simplify <- strats[strats[,2]>=0&strats[,1]>=0,]
strats_opt <- strats_simplify[strats_simplify[,4]==max(strats_simplify[,4]),]
TernaryPoints(c(strats_opt[2],strats_opt[3],strats_opt[1]),cex=1,col='red',bg='red',pch=21)
#strats_opt
Great, so now we have an optimization. Let’s put it into a function.
mixo_opt <- function(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI){
# Add analytical solution
strats <- matrix(data=NA,nrow=3,ncol=4)
## Strict phagotrophy
alpha_til_V <- max((qG*yGM+qC*r)/(qC*uCB*B + qG*yGM),(qG*yGM)/(qN*uNB*B+qG*yGM))
alpha_til_P <- 0
alpha_til_M <- 1 - alpha_til_V - alpha_til_P
g_til_CV <- qC*(uCB*B*alpha_til_V-r)
g_til_NV <- qN*uNB*B*alpha_til_V
g_til = qG*yGM*alpha_til_M # growth rate of strict phagotroph
strats[1,1] <- alpha_til_V
strats[1,2] <- alpha_til_P
strats[1,3] <- alpha_til_M
strats[1,4] <- g_til
## Strict phototrophy
alpha_bar_V = 0
alpha_bar_P = max(c((qG*yGM+qC*r)/(qC*uCL*L+qG*yGM),qG*yGM/(qN*uNI*I+qG*yGM)))
alpha_bar_M = 1-alpha_bar_P
g_bar_CP = qC*(uCL*L*alpha_bar_P - r)
g_bar_NP = qN*uNI*I*alpha_bar_P
g_bar = qG*yGM*alpha_bar_M # growth rate of strict phototroph
strats[2,1] <- alpha_bar_V
strats[2,2] <- alpha_bar_P
strats[2,3] <- alpha_bar_M
strats[2,4] <- g_bar
## Mixotrophy
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)
alpha_star_V <- (term.c-term.d)/(term.a+term.b)
alpha_star_P <- (qG*yGM*(1-alpha_star_V)+qC*r-qC*uCB*B*alpha_star_V)/(qG*yGM+qC*uCL*L)
alpha_star_M <- 1 - alpha_star_V - alpha_star_P
g_star = qG*yGM*alpha_star_M
strats[3,1] <- alpha_star_V
strats[3,2] <- alpha_star_P
strats[3,3] <- alpha_star_M
strats[3,4] <- g_star
## Picking optimum
strats_simplify <- strats[strats[,2]>=0&strats[,1]>=0,]
strats_opt <- strats_simplify[strats_simplify[,4]==max(strats_simplify[,4]),]
Res <- list("Optimum" = strats_opt)
}
# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,1))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L
# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01
# Uptakes
yGM <- 1 # yield from growth
uCL <- 0.01 # carbon from light
uNI <- 1 # nutrients from inorganic nutrients
uCB <- 0.000003 # carbon from bacteria
uNB <- uCB/10 # nitrogen from bacteria
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# External resource availability
L <- 10
I <- 2.8 #3 #2
B <- 1e6
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
# Compute the strategy for all investment strategies
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){ qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){ qG*yGM*b }
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){ values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }
# Add shading to the plot
#ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))
# Compute the growth-maximizing strategy
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
TernaryPoints(c(Res$Optimum[2],Res$Optimum[3],Res$Optimum[1]),cex=1,col='red',bg='red',pch=21)
Without dynamic resources, if the dilution rate is lower than the initial growth rate, the population will grow without bound.
# Parameters
# External resource availability
L0 <- 100
I0 <- 10
B0 <- 1e6
# Respiratory costs
r <- 0.1#0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Dilution rate of the chemostat
D <- 0.001
# Set up the simulation
tset <- seq(from = 0, to = 1000, by = 0.01)
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0
growth <- NaN*tset; growth[1] <- Res$Optimum[4]
# Run the simulation
for(i in 2:length(tset)){
dt <- tset[i] - tset[i-1]
L <- L0
B <- B0
I <- I0
jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
jG <- ( yGM*Mset[i-1] )
growth[i] <- min(qC*jC,qN*jN,qG*jG)
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
alphaV <- Res$Optimum[1]
alphaP <- Res$Optimum[2]
alphaM <- Res$Optimum[3]
Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaP - D*Pset[i-1])
Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaV - D*Vset[i-1])
Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaM - D*Mset[i-1])
}
Tset <- Pset+Vset+Mset
plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='')
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
Adding in chemostat dynamics for the mixotroph
# Parameters
# External resource availability
L0 <- 100
I0 <- 10
B0 <- 1e6
# Respiratory costs
r <- 0.1#0.01
# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Dilution rate of the chemostat
D <- 0.001
# Set up the simulation
tset <- seq(from = 0, to = 1000, by = 0.01)
Res <- mixo_opt(L0,I0,B0,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]*.5
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0
growth <- NaN*tset; growth[1] <- Res$Optimum[4]
# Run the simulation
for(i in 2:length(tset)){
dt <- tset[i] - tset[i-1]
L <- L0
B <- B0
I <- I0
jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
jG <- ( yGM*Mset[i-1] )
growth[i] <- min(qC*jC,qN*jN,qG*jG)
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
alphaV <- Res$Optimum[1]
alphaP <- Res$Optimum[2]
alphaM <- Res$Optimum[3]
Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaP - D*Pset[i-1])
Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaV - D*Vset[i-1])
Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaM - D*Mset[i-1])
}
Tset <- Pset+Vset+Mset
par(mar=c(1,4,3,2),mfrow=c(2,1))
plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='',xaxt='n')
axis(side=1,labels=F)
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
#mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
par(mar=c(4,4,0,2))
plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
And now dynamic resources
# Parameters
# External resource availability
L0 <- 100 #1000
I0 <- 10 #0#10
B0 <- 1e6
# Respiratory costs
r <- 1#0.1#0.01
# Uptakes
yGM <- 1
uCL <- 0.01 #0.01*2
uNI <- 1
attack_rate <- 0.000003
f <- 1 # conversion efficiency of bacterial C to mixotroph C
uCB <- attack_rate*f
uNB <- uCB/10 ##5*attack_rate*f
# Stoichiometries
qC <- 1
qN <- 1
qG <- 1
# Dilution rate of the chemostat
D <- 0.1
# Absorptivities of chemostat components
kappa0 <- 0 # base absorptivity of water
kappaI <- 0 # absorptivity of inorganic nutrients
kappaB <- 0.#00001 # per capita absorptivity of bacteria
kappaP <- 0.#001 # per capita absorptivity of plastids
kappaV <- 0.#00001 # per capita absorptivity of vacuoles
kappaM <- 0.#00001 # per capita absorptivity of growth machinery
# Set up the simulation
tset <- seq(from = 0, to = 20, by = 0.001)
Res <- mixo_opt(L0,I0,B0,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0
alphaVset <- Vset; alphaVset[1] <- Res$Optimum[1]
alphaPset <- Vset; alphaPset[1] <- Res$Optimum[2]
alphaMset <- Vset; alphaMset[1] <- Res$Optimum[3]
growth <- NaN*tset; growth[1] <- Res$Optimum[4]
# Run the simulation
for(i in 2:length(tset)){
dt <- tset[i] - tset[i-1]
# Determine resource acquisition
L <- Lset[i-1]
B <- Bset[i-1]
I <- Iset[i-1]
jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
jG <- ( yGM*Mset[i-1] )
# Compute growth rate
growth[i] <- min(qC*jC,qN*jN,qG*jG)
# Determine optimal investment strategy
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
alphaVset[i] <- Res$Optimum[1]
alphaPset[i] <- Res$Optimum[2]
alphaMset[i] <- Res$Optimum[3]
# Update mixotroph state variables
Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaPset[i] - D*Pset[i-1])
Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaVset[i] - D*Vset[i-1])
Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaMset[i] - D*Mset[i-1])
# Update resource state variables
Bset[i] <- Bset[i-1] + dt*( D*(B0 - Bset[i-1]) - attack_rate*Vset[i-1]*Bset[i-1] )
Iset[i] <- Iset[i-1] + dt*( D*(I0 - Iset[i-1]) - uNI*Iset[i-1]*Pset[i-1] )
Lset[i] <- L0*exp(-kappa0 - kappaP*Pset[i] - kappaV*Vset[i] - kappaM*Mset[i] - kappaB*Bset[i] - kappaI*Iset[i])
}
Tset <- Pset+Vset+Mset
par(mar=c(1,4,3,2),mfrow=c(2,1))
plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='',xaxt='n')
axis(side=1,labels=F)
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset, alphaPset, lty = 2, col = Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,alphaVset, lty = 2, col = Vcol)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
lines(tset,alphaMset,lty=2, col = Mcol)
#mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='right',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=F,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
par(mar=c(4,4,0,2))
plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='',ylim=c(min(Pset,Vset,Mset,Tset),max(Pset,Vset,Mset,Tset)),log='y')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
lines(tset,Tset,col='black',lwd=2)
#legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)
plot(tset,alphaPset,type='l')
plot(Bset,Iset,type='l',lwd=2)